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Abstract. Two types of surface models have been investigated by Monte Carlo simulations on triangulated 
spheres with compartmentalized domains. Both models are found to undergo a first-order collapsing transi- 
tion and a first-order surface fluctuation transition. The first model is a fluid surface one. The vertices can 
freely diffuse only inside the compartments, and they are prohibited from the free diffusion over the surface 
due to the domain boundaries. The second is a skeleton model. The surface shape of the skeleton model is 
maintained only by the domain boundaries, which are linear chains with rigid junctions. Therefore, we can 
conclude that the first-order transitions occur independent of whether the shape of surface is mechanically 
maintained by the skeleton (= the domain boundary) or by the surface itself. 
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1 Introduction 

The crumpling transition of membranes is an interesting 
topic in the softmatter physics as well as in the biological 
physics [1,2, 3]- A well-known model for such transition is 
the surface model of Helfrich, Polyakov, and Kleinert [H 
03G>] . A considerable number of theoretical and numerical 
studies have been devoted to reveal the phase structure 
of the model F7U51fWT01JTlirT^rT^rT^1JT31JT^rT^ . Recently, it 
was shown by Monte Carlo (MC) simulations that the 
model undergoes a first-order transition on spherical and 
fixed connectivity surfaces |llU12j . and the transition is 
universal [13j . The vertices can move only locally on the 
surface because of the fixed connectivity nature in those 
surface models. 

However, the crumpling transition is not yet clearly un- 
derstood in biological membranes. If we consider the pos- 
sibility of the transition in the cell membranes, we should 
take account of the fluid nature such as the lateral diffu- 
sion of lipids. 

Conventionally, the free diffusion of lipids has been 
realized by the dynamical triangulation technique in the 
surface model [14lfT5l[T6l[T7] . The diffusion of lipids has no 
cost in energy on fluid surfaces. 

In the cell membranes, however, the free diffusion of 
membrane proteins and lipids is suffered from heteroge- 
neous structures. Such molecules are known to undergo 
the so-called hop diffusion over the surface, which was re- 
cently found experimentally [18] . The free diffusion of the 
molecules is prohibited due to the cytoskeleton. The dif- 
fusion rate is, therefore, 10-100 times lower than that of 
artificial membranes, which are usually homogeneous and 
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have no such domain structure. Moreover, some artificial 
membranes are considered to have skeletons, because they 
are partly polymerized |19] . 

Motivated by this fact observed in the cell membranes, 
we study firstly in this paper a dynamically triangulated 
surface model with compartmentalized domains whose bound- 
aries are composed of triangle edges (or bonds) that are 
not to be flipped. The diffusion is constrained so that ver- 
tices can diffuse only inside the compartments, and hence 
the vertices never jump across from one compartment to 
the other compartments. Nevertheless, we consider that 
such constrained lateral diffusion can simulate the hop dif- 
fusion in the cell membranes as the first approximation. 

It is also interesting to see whether the transition oc- 
curs in a simplified skeleton model, where only skeletons 
maintain the mechanical strength of the surface. Skeleton 
models for the cytoskeleton were already investigated in 
[2"0] . A hard- wall and hard-core potential was assumed on 
the polymer chains with junctions, and the responses to 
some external stress and the compression modulus were 
obtained 20 . Giant fluid vesicles coated with skeletons 
was experimentally investigated, and the mechanical prop- 
erties were reported |21) . where the actin filaments intro- 
duce an inhomogeneous structure in the homogeneous ar- 
tificial membrane. 

In [52], the phase structure of a surface model with 
skeleton was investigated, and it was reported that the 
model has a first-order transition between the smooth 
phase and the crumpled phase. The interaction of the 
model in |22j is described by a one-dimensional bending 
energy for linear chains (or bonds) and the two-dimensional 
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bending energy for junctions. The two-dimensional Gaus- 
sian bond potential is also assumed in the Hamiltonian. 
Consequently, a simplified skeleton model can be obtained 
from the model in [22] by replacing the elastic junction 
with a rigid junction. 

Therefore, we study secondly in this paper the rigid 
junction model by using the canonical Monte Carlo sim- 
ulations and see how the transition of [22j occurs in such 
a simplified model. We must note that the rigid junction 
model is not identical to the elastic junction model in [22 . 
In fact, there are two types of elasticity at the junctions; 
one is the out-of-plane elasticity and the other is the in- 
plane elasticity. The former elasticity can be rigid in the 
limit of infinite bending rigidity bj — > oo in the elastic junc- 
tion model of 22J, however, the in- plane elasticity can not 
be controlled in the elastic junction model. Therefore, the 
rigid junction model in this paper and the elastic junction 
model in [22] are considered to be two different models. 



2 Fluid Surface Model 

By dividing every edge of the icosahedron into £ pieces 
of the uniform length, we have a triangulated surface of 
size N = 10£ 2 + 2 (= the total number of vertices). The 
starting configurations are thus characterized by iV 5 = 12 
and Nq — N— 12, where N q is the total number of vertices 
with the co-ordination number q. 

The compartmentalized structure is built on the sur- 
face by keeping the boundary bonds unflipped in the MC 
simulations with dynamical triangulation. The boundary 
of the compartment is constructed from a sequence of 
bonds that remain unflipped. The total number Nq of 
compartments depends on the surface size N. We fix n 
the total number of vertices inside a compartment to the 
following values: 

n = 21, 36, 66, 91, 120, (1) 
(# of vertices in a compartment). 

As a consequence, Nc is increased with the increasing N. 
The reason why we fix n is that the size of compartment is 
considered to be finite, and then it is expected that total 
number of lipids in the compartment also remains finite in 
the cell membranes. We must emphasize that the finitcness 
of n, rather than the value of n, is physically meaningful, 
because we do not always have one to one correspondence 
between the vertices and the lipid molecules. 

Figures [Ha), (b) show surfaces of {N,n) — (2562,21) 
and {N,n) = (15212,66) for the starting configurations 
of MC simulations. Thick lines denote the compartment 
boundary consisting of the bonds that are not to be flipped. 
Vertices on the boundary of compartments can locally 
fluctuate, and they are prohibited from the diffusion. The 
remaining vertices freely diffuse only inside the compart- 
ments. The starting configurations in Figs.[lja),(b) for the 
fluid model simulations are almost identical to those for 
the skeleton model simulations, which will be defined in 
the following section. 




(a) (N, n) = (2562, 21) (b) (TV, n) — (15212, 66) 



Fig. 1. Starting configuration of surfaces of (a) (N,n) = 
(2562,21) and (b) (N,n) = (15212,66), where n is the total 
number of vertices inside one compartment. Thick lines denote 
the compartment boundary consisting of bonds. 

We note that the fixed connectivity surface model is 
obtained from the compartmentalized fluid model in the 
limit n — ► 1, where the vertices are prohibited from the 
free diffusion. On the contrary, we obtain the fluid surface 
model in the limit of N , where all the vertices freely 
diffuse over the surface. 

The compartmentalized fluid surface model is defined 
by the partition function 

' ,/ n 

T J i=l 

S(X,T) = S 1 +bS 2 , 

where b is the bending rigidity, J denotes that the center 
of the surface is fixed in the integration. S(X, T) denotes 
that the Hamiltonian S depends on the position variables 
X of the vertices and the triangulation T. Y^t denotes 
the sum over all possible triangulations T, which keep the 
compartments unflipped. The Gaussian term Si and the 
bending energy term S*2 are defined by 

S 1 =^2(X i -X j )\ S 2 = - m • n,), (3) 

(y) 

where Y(ij) m $i is the sum over bonds (ij) connecting 
the vertices i and j, and Y^Uj) m ^2 is also the sum over 
bonds (ij), which are edges of the triangles i and j. The 
symbol rij in S% denotes the unit normal vector of the 
triangle i. We emphasize that the compartment boundary 
gives no mechanical strength to the surface in this model. 

The bending rigidity b has unit of kT, where k is the 
Boltzmann constant, and T is the temperature. The sur- 
face tension coefficient a of Si is fixed to be a = 1; this 
is always possible because of the scale invariant property 
of the model. In fact, in the expression aSi + bS2 we im- 
mediately understand that a = 1 is possible, because the 
coefficient a of Si can be eliminated due to the scale in- 
variance of the partition function. Since the unit of a is 
(1/length) 2 , the length unit of the model is given by yT/a. 
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We use the unit of length provided by yl/a = 1 in this 
paper, although a is arbitrarily chosen to be fixed. 



3 Skeleton Model 

The model is defined on a triangulated surface, which is 
characterized by N the total number of vertices includ- 
ing the junctions, N$ the total number of vertices on the 
chains, Nj the total number of junctions, and L the length 
of chains between junctions. The junctions are assumed as 
rigid plates; twelve of them are pentagon and the others 
are hexagon. It should be noted again that Nj is included 
in TV; a junction is counted as a vertex. 

The surface of size (N,N s ,Nj,L) = (2322,600,42,6) 
corresponds to that shown in FigQJa) for the fluid model. 
The reason why N — 2322 of the surface is smaller than 
N = 2562 of the one in FigQTa) is because the junctions are 
rigid objects. One hexagonal junction reduces N by 6, and 
one pentagonal junction also reduces iV by 5 if they were 
assumed as rigid objects. Thus, we can check that 2562 = 
2322+6 x 30+5 x 12, where 30 and 12 are the total number 
of pentagonal junctions and that of hexagonal junctions, 
respectively. 

We fix the chain length L such that 

L = 6(n=21), L = 11 (ti = 66), (4) 

which respectively correspond to the values n = 21, n = 
66; the total number of vertices inside a compartment in 
Eq.JT]). The reason why we fix L is the same as that for n. 

The Hamiltonian of the model is given by a linear com- 
bination of the two-dimensional Gaussian bond potential 
S\ and the one-dimensional bending energy S% > which 
are defined by 

St^iXi-Xtf, =J^[l-cose m ]. (5) 

(y) to) 

In these expressions, Ylnj) m denotes the sum over all 
the bonds (ij) connecting the vertices i and j, and X+y) 

in S% denotes the sum over bonds i and j, which contain 
not only bonds in the chains but also virtual bonds that 
connect the center and the corners of the rigid junctions. 
The symbol 6uj-\ in is the angle between the bonds 
i and j, which include the virtual bonds. The Gaussian 
potential Si is defined on all the bonds including those on 
the chains. As a consequence, the model is considered to 
be a surface model, although the mechanical strength is 
maintained by one-dimensional clastic skeletons joined to 
each other at the rigid junctions. 

Triangulated spherical surfaces for the skeleton model 
are almost identical to those shown in FigsH](a),(b) as 
mentioned above. Only difference between the surfaces 
is whether the junctions are rigid objects or not. Figure 
[2] shows a hexagonal rigid junction connected to chains, 
where the angle is defined not only at the vertices 
on the chains but also at the corners (=virtual vertices) 




Fig. 2. A hexagonal junction connected to chains. The angle 
9uj) in 5*2 is defined not only at the vertices on the chains 
but also at the corners (=virtual vertices) of the junction. The 
triangular lattices attached to the chains were eliminated from 
the figure to clarify the chains. 



of the junction. The triangular lattices attached to the 
chains were eliminated from Fig. [2] to clarify the chains. 

We should comment on the size of the junctions. The 
junctions are two-dimensional objects and therefore have 
their own size to be fixed. The size of junction can be 
specified by the edge length R; the perimeter length of the 
pentagonal (hexagonal) junction is therefore expressed by 
5R (6R). In this paper, we fix the size of the junctions 
such that 

R = 0.1. (6) 

The value R = Q.l is quite smaller than that of the elas- 
tic junctions in [22], where the edge length squared is 
R 2 ~ 0.5 because of the relation Si/N = 1.5, which is 
also satisfied in the fluid model defined in the previous 
section. 

We must note that the junction size in Fig. [2] were 
drawn larger, comparing to the bond length, than that 
expected from Eq.©. In the following section, we discuss 
how do we fix the size R to the assumed value in Eq.® 
in the MC simulations. 

The partition function Z of the skeleton model is de- 
fined by 

,/ N 

Z= Y[dX t exp[-S(X)}, S(X) = Si + bS^\ (7) 

i— 1 

where b is the bending rigidity corresponding to the one- 
dimensional bending energy, and J denotes that the cen- 
ter of the surface is fixed. The integration YiiLi dXi is a 
product of the integration over vertices and that of junc- 
tions such that 

N 

Y[dXi= JJ dX t JJ dX l7 (8) 

i— 1 vertices i junctions i 

where functions i ^ s t ne integration over the degrees 
of freedom for the three-dimensional translations and ro- 
tations. 
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4 Monte Carlo technique 

The vertices X are shifted so that X' — X+SX , where SX 
is randomly chosen in a small sphere. The new position X' 
is accepted with the probability Min[l, exp(— AS)], where 
AS = S (new) -S (old). 

The summation over T in the fluid model partition 
function Z of Eq.((2]) is performed by the standard bond 
flip technique. The bonds are labeled with sequential num- 
bers. The total number of bonds Nb is given by Nb — 
37V— 6, which includes the bonds in the boundary of com- 
partments. Firstly, the odd-numbered bonds are sequen- 
tially chosen, to be flipped and secondly, the remaining 
even-numbered bonds are chosen. The flip is accepted with 
the probability Min[l, exp(— AS)]. In this procedure, the 
compartment boundary remains unfiipped. N updates of 
X and Nb/2 updates of T are consecutively performed 
and make one MCS (Monte Carlo Sweep). The radius of 
the small sphere for SX is chosen so that the rate of accep- 
tance for X is about 50%, which is controlled by a small 
number for the radius SX at the beginning of the simu- 
lations. We introduce the lower bound 1 x 1CP 8 for the 
area of triangles. No lower bound is imposed on the bond 
length. 

The assumed sizes in the fluid model simulations are 
listed in Table [T] Three sizes of surfaces are assumed for 
each compartment size n except for n = 120. The third 
size for n — 91 is relatively large and therefore time con- 
suming for the fluid simulations, and the size N = 11562 is 
sufficiently large to show that there is no phase transition 
on the surface of n = 120. 

Table 1. The size N of surfaces for the fluid model simulations. 
Three sizes of surfaces are assumed for each compartment size 
n except for n= 120. 

n = 21 w = 36 w = 66 n = 91 n = 120 

2562 1002 1692 2252 2892 

5762 4002 6762 9002 11562 

10242 9002 15212 20252 



Total number of MCS after the thermalization MCS 
at b close to the transition point is about 2xl0 8 ~ 3xl0 8 
MCS on the N > 9002 surfaces and 1 x 10 8 - 1.7 x 10 s 
MCS on the N < 6762 surfaces, and relatively smaller 
number of MCS at b far from the transtion point. The 
thermalization MCS 1 x 10 7 - 1.7 x 10 8 MCS on the iV> 
9002 surfaces and 1 x 10 7 - 3 x 10 7 MCS on the N < 6762 
surfaces. The reason of such a large number (1.7 xlO 8 ) 
of thermalization MCS seems due to the discontinuous 
nature of the transition. The starting configurations of MC 
are just like those shown in FigsUla) andQJb), and hence 
they are in the smooth phase. In fact, large surfaces in 
the collapsed phase close to the transition point eventually 
collapsed after such a long thermalization MCS. 

The update of X in MC for the skeleton model par- 
tition function in Eq.Q can be divided into two steps, 



which are corresponding to the integrations Ilvcrticos t ^ i 
an d rij unctions j dXi in Eq.®. The first is the update of 
X of vertices including those in the chains. The second 
is the update of the position of the junctions as three- 
dimensional rigid objects. This can be further divided 
into two processes: the first is a random three-dimensional 
translation, and the second is a random three-dimensional 
rotation. All of these MC processes are independently per- 
formed under about 50% acceptance rate. 

The junction size R is fixed to i? = 0.1 in Eq.Q during 
the thermalization MCS. The initial value of R is given by 
R ~ 0.7 on the surfaces such as those shown in FigsQTa) 
and QTb). Thus, we reduce R from R ~ 0.7 to R = 0.1 
by AR = 6 x 10~ 6 at every 25 MCS in the first 2.5 x 10 6 
MCS. Because of this forced reduction of the junction size, 
the equilibrium statistical mechanical condition seems to 
be violated in the first 2.5 x 10 6 MCS. Therefore, rela- 
tively many thermalization (1.75 x 10 7 or more) MCS is 
performed after the first 2.5 x 10 6 MCS for the reduction. 

We use surfaces of size listed in Table [5] for the skeleton 
model simulations. Four sized are assumed for L = 6, and 
three sizes for L = ll. 



Table 2. The size N of surfaces for the skeleton model simu- 
lations. Four sizes of surfaces are assumed for L — 6, and three 
sizes for L = ll. 
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6522 
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6 


9282 


2400 
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11 


14672 


2700 


92 


6 


14502 


3750 


252 


11 


26082 


4800 


162 


6 


20882 


5400 


362 











The total number of MCS is about 1.6 x 10 8 ~ 2x 10 8 
for the N < 6522 surfaces and 3 x 10 8 ~ 4 x 10 8 for the 
N > 9282 surfaces. The thermalization MCS is 2 x 10 7 for 
the N < 6522 surfaces and 2xl0 7 - 1x10 s for the N < 9282 
surfaces. 

A random number sequence called Mersenne Twister 
1231 is used in the simulations. 



5 Fluid Model Simulations 

First, we show in Figs. EHa),(b) snapshots of the (N,n) = 
(20252,91) surfaces obtained at 6= 1.57 in the collapsed 
phase and at b— 1.58 in the smooth phase. Figures[3^c),(d) 
show the surface sections. Thus, we confirmed that the 
smooth phase can be seen at finite b close to the transition 
point. 

The crumpling transition is conventionally understood 
as the one of surface fluctuation accompanied by surface 
collapsing phenomena, which can be seen in our surface 
model as we have seen in the snapshots in Figsi3][a)-[3jd). 
Therefore, we expect that the mean square size X 2 reflects 
the collapsing transition on spherical surfaces. The mean 



(a) (20252,91), 6=1.57 
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(b) (20252,91), 6=1.58 



value corresponding to the transition point of the fixed 
connectivity surface model p~2] • We can also confirm that 
b n disappears in the limit of n — > N at sufficiently large 
N, because we find no transition at n = 120. 




(c) The surface section 



(d) The surface section 



Fig. 3. Snapshot of the fluid surface of (N, n) = (20252, 91) 
obtained at (a) 6 = 1.57 in the collapsed phase and at (b) 
6=1.58 in the smooth phase, and (c), (d) the surface sections. 
The mean square size X 2 defined in Eq.® is X 2 ~ 15 in (a) 
and X 2 — simeql5l in (b). 



1.6 1.65 

(a) " (b) " (c) " (d) " (e) b 

Fig. 5. The bending energy S2/NB versus 6 obtained on the 
surface of (a) n = 21, (b) 71 = 36, (c) n = 66, (d) 71 = 91, and (e) 
n = 120. Nb(^3N — 6) is the total number of bonds. 



The crumpling transition is originally understood as 
the one for surface fluctuation phenomena. Therefore, the 
bending energy S 2 /N B , defined in Eq.©, is expected to 
reflect the transition. Figures Efa),(b),(c),(d) and (e) show 
S%/Nb versus b corresponding to n — 21, n — 36, 71 = 66, 
n = 91, and n — 120, respectively. Discontinuous change 
of S2/NB can also be seen at b where X 2 discontinuously 
changes, although it is not sufficiently clear in the figures. 
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(a) " (b) " (c) " (d) 

Fig. 4. The mean square size X 2 versus 6 obtained on the fluid 
surfaces of (a) 71 = 21, (b) n = 36, (c) n = 66, (d) 71 = 91, and 
(e) n — 120. Solid lines were obtained by the multihistogram 
reweighting technique. 



square size X 2 is defined by 



X 2 



- y & 



X) 



X 



N ^ 



(9) 



where X is the center of the surface. FiguresHfa),(b),(c),(d), 
and (e) show X 2 obtained on the surfaces of n = 2l, n = 36, 
n = 66, n = 91, and n = 120, respectively. Solid lines 
connecting the data were obtained by the multihistogram 
reweighting technique [24J . A discontinuous change of X 2 
can be seen in the cases of n < 91 when the size N in- 
creases, while X 2 largely fluctuates at n — 120. This in- 
dicates that a first-order transition occurs at n < 91 and 
that it disappears at n = 120. We also find that the tran- 
sition point b n moves left on the 6-axis as n decreases. It 
is expected in the limit of n — * 1 that b n reduces to the 




Fig. 6. The variation of the gap of S2/NB against n, which 
were obtained on the surfaces of size (N,n) = (10242,21), 
(9002, 36), (15212, 66), (20252, 91), and (11562, 120). The solid 
circles (•) and the dashed lines denote the value of Sz/Nb in 
the smooth (the crumpled) phase at the first-order transition 
point. The diamonds (©) correspond to the results obtained on 
the surface (20252,91). The symbol (®) denotes the critical 
point of transition, where the discontinuous transition termi- 
nates and turns to a continuous or a higher-order one. The 
critical value n c is expected to be n c — 91. 



In order to see the gap of S2/NB more clearly, we show 
the variation of the gap of S2/NB against n in FigJB] The 
solid circles (•) and the dashed lines denote the value of 
S2/NB in the smooth (the crumpled) phase at the first- 
order transition point. The diamonds (o) correspond to 
the results obtained on the surface (20252,91). We find 
that the value of S2/NB in the smooth phase (or in the 
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crumpled phase) increases as n decreases at the transi- 
tion point. The value of S2/NB at the transition point 
becomes identical to the one of the fixed connectivity sur- 
face model in the limit of n — > 1 [12] . On the contrary, the 
discontinuous change of S^/Nb is expected to disappear 
at sufficiently large n, because the discontinuity of S2/NB 
seen clearly at n = 21 ~ 66 disappears when n— >N at suf- 
ficiently large N . This implies that there exists a finite n c , 
where the first-order transition terminates and turns to 
a continuous or a higher-order one. The gap of S2/NB at 
n = 21 reduces as n increases and eventually goes to zero at 
n = n c , which is expected to be n c ~91. At n=91, S2/NB 
seems continuous while X 2 is discontinuous as confirmed 
in Figged). 

We note that the maximum co-ordination number q max 
obtained throughout the simulations is as follows: q m & x = 
38 on the (N, n) = (10242, 21) surface at 6 = 1.38, g max = 42 
on the (iV,n) = (9002,36) surface at 6=1.42, g max = 40 on 
the (N, n) = (15212, 66) surface at 6 = 1.51, q max = 42 on 
the (N, n) = (20252, 91) surface at 6=1.55, and g raax = 44 
on the (N,ri) = (11562, 120) surface at 6=1.56. 



because of the fact that the surface of (20252, 91) becomes 
smooth at 6 > 1.58 and crumpled at 6 < 1.57, which was 
clarified from X 2 in Fig[4jd). Thus, the critical value was 
expected to be n c ~91. 




Fig. 8. (a) X 2 vs. b, (b) S2/NB versus b, and (c) Cs 2 versus 
b obtained on the fluid surface model without the compart- 
mentalized structure. The symbol n = N drawn on the figures 
denotes that the surfaces are those without the compartment. 
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Fig. 7. The histogram h(S<2) and the corresponding vari- 
ation of S2/NB against MCS obtained on the surfaces of 
(a),(b) (N,n) = (5762,21), (c),(d) (N, n) = (4002, 36), (e),(f) 
(JV,n) = (6762,66), and (g),(h) (N, n) = (9002, 91). 



In order to show the discontinuity in S2/NB more 
clearly, we plot in Figs[7]Ja)-[7fh) the distribution (or his- 
togram) h{S2) of S2/NB and the corresponding varia- 
tion of S2/Nb- These were obtained on the surfaces of 
(N,ri) = (5762,21), (4002,36), (6762,66), and (9002,91). 
The discontinuity of S^/iVs can be seen in the histogram 
on the N > 5762 surfaces. Because of the size effect, the 
transition appears to be continuous on the N = 4002 
surface. The double peak at (N, n) — (5762, 21) is more 
clear than that at (N,n) — (6762,66), because the gap 
of S2/NB reduces as n increases. We should note that 
the double peak structure is very hard to see in h(S2) on 
the N > 9002 surfaces at n = 36. When the configura- 
tion is once trapped in the smooth (the crumpled) state, 
it hardly changes to the crumpled (the smooth) state at 
the transition point on such large surfaces. This problem 
may be resolved with more sophisticated MC techniques 
[27]l28] . No double-peak structure is found in Fig[7]^g) and 
the variation of S2/NB smoothly varies in FigJTJh). This 
indicates that the transition is a continuous one at n=91 



The phase transition is expected to disappear from the 
fluid surface model defined on the surfaces that have no 
compartment. In order to show this, we performed MC 
simulations on the surfaces without the compartments up 
to the size TV = 4842. Figure Ufa) shows X 2 versus 6. We 
can see no abrupt growing of X 2 in the figure. The bending 
energy S2 /Ng shown in Fig. [S^b) , where the variation of 
52 versus 6 seems almost independent of N. The specific 
heat Cs 2 , which is defined by C S2 = {b 2 /N)( (5 2 -(S" 2 )) 2 ), 
is expected to reflect the phase transition. However, we 
can see no anomalous behavior in Cs 2 shown in Fig. [Hf c) ; 
there can be seen no peak in Cs 2 ■ Thus we confirmed that 
the phase transition disappears from the model if n—*-N 
at sufficiently large N. 

We must comment on the relation between the above 
results and those of fluid surface simulations in [25], be- 
cause the phase transition can be seen in [25j while it can 
not in Figs[SJa),(b). The difference between the bond flip 
procedure in this paper and that of [25] is the reason why 
the phase transition can be seen in the model in [25] and 
it can not be seen in the same model in Figs[H]Ja),(b). In 
the simulations of this paper we label the bonds by a se- 
quence of numbers as stated above and perform the bond 
flip by using the sequential numbers. On the contrary, a 
vertex is firstly chosen randomly in [25], and secondly a 
bond is randomly chosen to be flipped from the bonds 
connecting the chosen vertex. As a consequence, this pro- 
cedure gives a large (small) weight to the vertices which 
have small (large) co-ordination number in the dynami- 
cal triangulations. Therefore, the procedure in [25j seems 
change effectively the coefficient a of the co-ordination 
dependent term — cvj^ i logc^, which comes from the inte- 
gration measure Yij dXjqf, where a is fixed to a = in 
this paper and in [35]. We know that the phase structure 
depends on the co-ordination dependent term in the fluid 
surface model 1261. 
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6 Skeleton Model Simulations 

Snapshots of the skeleton surfaces are shown in FigsGSa)- 
EKd). Figure Efa) is a surface of size (N,N s ,Nj,L) = 
(26072,4800,162,11) obtained at b = 12.3 in the crum- 
pled phase, and FigO^b) is the one obtained at b— 12 A in 
the smooth phase. The surface sections of FigsJUa) and 
EJb) are shown in FigsEJc) and[9jd), respectively. These 
figures were drawn in the same scale. We immediately see 
the surface in the smooth phase is actually swollen, while 
the surface is collapsed in the collapsed phase. 



(a) 6=12.3 (Collapsed) 



(b) 6=12.4 (Smooth) 




(c) The surface section 



(d) The surface section 



Fig. 9. Snapshot of the surface of size (N,Ns,Nj,L) — 
(26082,4800,162,11) obtained in the collapsed phase at (a) 
6= 12.3 and in the smooth phase at (b) 6= 12.4, both of which 
are close to the transition point. The mean square size X 2 is 
X 2 ~12 in (a) and X 2 ~98 in (b). 



The Gaussian bond potential Si/N is shown against b 
in Figs fTUF a) and fTOT b). which correspond to the lengths 
L = 6 and L = ll, respectively. The values of Si/N in the 
figures slightly deviate from Si/N = 1.5, which is satis- 
fied on the surface without the rigid junctions or the rigid 
junctions of negligible size. The reason of this discrepancy 
is because the surface includes the rigid junctions of fi- 
nite size. A vertex is the zero-dimensional point, while the 
rigid junction is a two-dimensional plate and hence shares 
an area of the surface. We find a gap or a jump in Si/N 
of the (N,N s ,Nj,L) = (26082,4800,162,11) surface in 
FigfTOFb), which can be viewed as a sign of a discontinu- 
ous transition. 

Figures [TTTa) and flTTb) are plots of X 2 against b ob- 
tained under L = 6 and L = ll. We find that the variation 
of X becomes rapid against b as N increases. Thus, it is 



S,/N 



1.518 



1.514 



L=6 



O :N=20882 
A :N=9282 
□ :N=5222 



Sj/N 
1.512 

1.508 



L=ll 



O :N=26082 
A:N=14672 
□ :N=6522 



5.2 , 5.4 

b (b) 



12.5 , 13 
b 



(a) 

Fig. 10. The Gaussian bond potential Si/N versus 6 obtained 
on the surfaces of (a) L = 6 and (b) L = 11. Si/N slightly 
deviates from Si /N ~ 1.5. 
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N=5222 ■ 


5 


5.2 


b 54 




(a) 

Fig. 11. The mean square size X 2 versus 6 obtained on the 
surfaces of (a) L = 6 and (b) L=ll. The curves are drawn by 
the multihistogram reweighting technique. 



expected that the variation of X 2 has a jump at interme- 
diate value of b in either case of L. 

The bending energy /N' s is expected to reflect the 
transition, where N' s is given by 



N' s 



6Nj - 12. 



(10) 



because N' s 



S2 /N' s is the bending energy per vertex, 

is the total number of vertices where S^ 1 ' is defined, 
includes virtual vertices which are the corners of the junc- 
tions (see also Figj2]), which are not counted as vertices 
and hence are not included in JV5. Total number of vir- 
tual vertices are 6Nj — 12, because the hexagonal junction 
contains 6 virtual vertices, and the total number of pen- 
tagonal junction is 12. Thus, we have Eq. pU)) for N' s , and 
therefore we have N' s = 1890, N' s = 3360, N' s = 5250, 
and N' s — 7560 for the length L = 6 surfaces of size 
(N, N s , Nj) = (5222, 1350, 92), (9282, 2400, 162), (14502, 3750, 252), 
and (20882,5400,362); and N' s = 1440, N' s = 3240, and 
N' s = 5760 for the length L = 11 surfaces of size (N, Ns, Nj) = 
(6522,1200,42), (14672, 2700, 92), and (26082,4800,162). 

Figures fT2T a) and [T27 b) are plots of S2 /N' s against 
b obtained under L = 6 and L = ll. We find the expected 

behavior in /N' s under both conditions L — 6 and L — 

11; S^/N'g has a gap (or a jump) at intermediate b. This 
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clearly shows that the surface fluctuation transition is of 
first order. 



The two-dimensional extrinsic curvature S2 /N_ 



b is 



S 2 (1) /N s ' 



S 2 (1) /N s ' 



0.25 



0.23 



0.21 





1 1 1 1 1 




L=6 | ■ 




- N=5222 




V N=9282 


- X 




. N=20882 





0.11 



5.2 



5.4 



(a) 

Fig. 12. The one-dimensional bending energy S2 /N'§ versus 
6 obtained on the surfaces of (a) L = 6 and (b) L = ll. 




plotted in Figs fT5T a) and [T5f b) against b, which are cor- 
responding to the lengths L — 6 and L = ll. We find that 
the dependence of Sv; /Nb on b shown in Fig[T3]is almost 

identical to that of /N' s in FigfT2l The gap (or jump) 

(2) 

seen in Sv; /Nb also supports that the surface fluctuation 
transition is of first order. 

The specific heat corresponding to the one-dimensional 
bending energy is defined by 



0.09 



N: 



(12) 



which can also reflect phase transitions if it has an anoma- 
lous behavior. Figures fT4T a) and fT4T b) show C„(i) versus 

2 

b obtained under L — 6 and L=ll. Solid curves drawn in 
the figures were obtained by the multihistogram reweight- 
ing technique, and those curves apparently show an ex- 
pected anomalous behavior indicating that C„(i> is diver- 
The transition can also be reflected in the two-dimensiona^ en ^- wnen jy' qq ( or equivalently N — » 00) 
extrinsic curvature, which is defined by 



(ij) 



1 - ILi ■ Uj) , 



(11) 



where is the unit normal vector of the triangle i. The 

:(2) 



definition Eq. (fTT]l of S^ 2 ' is identical to that of S2 in 



Eq.f3|). In S 2 , Ysiij) denotes the summation over all 
nearest neighbor triangles i and j that have the common 
bond (ij), which includes bonds belonging to the skele- 
ton chains. Nb denotes the total number of bonds where 
S2 is defined, and it is given by NB = J2(ij) 1- Note that 

Nb includes virtual bonds, which are the edges of the rigid 

(2) 

junctions. In fact, we define S 2 even on the virtual bonds, 
because it is reasonable to define extrinsic curvature on 
those edges. It is also noted that S 2 is not included in 

(2) 

the Hamiltonian, and therefore Sv; gives no mechanical 
strength to the surface. 



'S2(l) 











1 | L=6 


N=20882 
,N=9282 - 

A 




o\ 

N=5222 



5.2 



5.4 



(a) b 

Fig. 14. The specific heat C 



-S2(l) 



40 



• N=26082 

N=14672 

j N=6522 




versus b obtained on 



the surfaces of (a) L — 6 and (b) L — 11. C (i) is defined by 



Eq. (|12p . The error bars on the symbols are the statistical error, 
which is obtained by the binning analysis. 




Fig. 13. The two-dimensional extrinsic curvature Sv; /Nb 
against b obtained on the surfaces of (a) L = 6 and (b) L = ll. 

is defined by Eq. Hllj) and is not included in the Hamilto- 

(2) 

nian. Nb is the total number of bonds where Si, is defined. 



The specific heat corresponding to the extrinsic curva- 

f 21 

ture S 2 > in Eq.JITj can also be defined by 

C s (2,=l<(^ 2) -<4 2) }) 2 ), (13) 
which reflects the transition as C„(i) does. Curvature co- 

(2) 

efficient for C„<2) was assumed to be 1, because S 2 is 

not included in the Hamiltonian. Figures [TBT a) and [TBT b) 
show C„<2) against b obtained under L — 6 and L = ll. We 

can see in C„(2j the same anomalous behavior as in C„(i) . 

°2 *2 

In order to see the anomalous behaviors in C„m and 

2 

C c (2) in more detail, we plot the peak values of them in 
2 

Figs fTBT a) and [TBTb) in log-log scales against N' s and N, 
respectively. The straight lines were drawn by fitting the 



H.Koibuchi: Phase transition of compartmentalized surface models 



9 



Cs2(2) C S2 (2) 




(21 

Fig. 15. The specific heat C (2) for S 2 versus b obtained on 

the surfaces of (a) L — 6 and (b) L — 11. C (2) is defined by 

Eg . (| 13[) . The error bars on the symbols are the statistical error, 
which is obtained also by the binning analysis. 



Cmax max 
S2(l) <-S2(2) 




against TV' obtained on the surfaces of L = 6 and L — 8. The 
straight lines are drawn by fitting the largest three data of 
and C™ to Eq.dHJ. The peak values and the statisti- 



cs a 2 

cal errors for the fittings were obtained by the multihistogram 
reweighting. 



data to 

c^w 1 - c^oc(ivr, (w) 

°2 °2 

where a\, 01 are critical exponents. Largest three data 
were contained in the fitting in the case L = 6. Thus, we 
have 

cr x = 0.911 ± 0.118, cr 2 = 0.608 ± 0.076, (L = 6), 
(7i = 1.51 ± 0.08, cr 2 = 0.948 ± 0.262, {L = 11). (15) 

The result 02 = 0.608(76) for L — 6 is inconsistent to the 
fact that the surface fluctuation transition is of first-order, 
however, <j\ = 0.911(118) is consistent to the discontinu- 
ous collapsing transition. The results u\ — 1.51(0.08) and 
a-2 = 0.948(262) under L = 11 support the discontinuous 
transition of surface fluctuation. 



7 Summary and Conclusion 

We have studied a compartmentalized fluid surface model 
and found that the model undergoes a first-order collaps- 
ing transition and a first-order surface fluctuation transi- 
tion between the smooth phase and the collapsed phase. 
The model is classified as a fluid surface model, although 
the long-range order and the phase transition can be seen 
at finite b. The compartmentalized structure is considered 
to be a reason for the existence of the phase transition. 
Consequently, the result is not in contradiction with the 
standard argument for the non-existence of long-range or- 
der in fluid membranes. Moreover, the critical point of the 
transition is strongly expected at finite n c (~ 91), where 
the collapsing transition and the surface fluctuation tran- 
sition terminate and turn into continuous or higher-order 
transitions. In fact, we demonstrated that the collapsing 
transition remains first-order at n < 91 and disappear at 
n= 120 on sufficiently large surfaces. It was also demon- 
strated that the phase transition of surface fluctuation re- 
mains first-order at n < 66 and disappear at n = 120 on 
sufficiently large surfaces, and that the transition weakens 
and still survives at n = 91. 

A surface model with skeletons has also been inves- 
tigated by using the canonical Monte Carlo simulations. 
The skeletons are composed of one-dimensional linear chains 
and rigid junctions, whose size is chosen sufficiently small 
compared to the mean bond length. The surface is a trian- 
gulated sphere and divided into a lot of compartmental- 
ized domains, whose boundary corresponds to the skele- 
tons, and it is almost identical to that for the compartmen- 
talized fluid model. The skeleton surface is characterized 
by (N,Ns,Nj,L), which are respectively the total num- 
ber of vertices including the junctions, the total number 
of vertices on the chains, the total number of junctions, 
and the length of chains between junctions. The length of 
chains was fixed to L = 6 and L = 11, which correspond 
to n = 21 and n = 66 the total number of vertices in a 
compartment. 

The mechanical strength is given to the surface only 
by the skeletons. There is no two-dimensional curvature 
energy in the Hamiltonian, while one-dimensional bending 
energy is defined on the compartment boundary. The two- 
dimensional Gaussian bond potential is included in the 
Hamiltonian just like in the standard surface model of 
Hclfrich, Polyakov and Kleinert. The skeleton model in 
this paper is different from the one with elastic junctions 
in [22], because the rigid junctions cannot be identified 
with the elastic junctions due to the property on the in- 
plane elasticity at the junctions. 

We found that the skeleton surface in this paper un- 
dergoes a first-order collapsing transition and a first-order 
surface fluctuation transition between the smooth phase 
and the crumpled phase. The one-dimensional bending en- 
ergy has a gap (or a jump) at intermediate bend- 
ing rigidity b, and the two-dimensional extrinsic curva- 
(2) 

ture S2 also has a gap at that point. These imply that 
the surface fluctuations are considered to be a first-order 
transition. Moreover, it is found that the mean square size 
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X 2 also has a gap at the transition point. This implies 
that the surface-collapsing phenomenon can be viewed as 
a first-order transition. 

The results in this paper together with those in [22] 
show that the first-order transitions can be seen in the 
spherical surface model even when the mechanical strength 
is maintained only by skeletons, which are composed of 
linear chains joined to each other at the junctions. More- 
over, the order of transitions is independent of whether 
the junction is elastic or rigid. 

We have studied two types of compartmentalized sur- 
face models; the mechanical strength is maintained by the 
two-dimensional curvature in the first model (the fluid 
model), and it is maintained by the one-dimensional cur- 
vature in the second model (the skeleton model). There- 
fore, we can also conclude that the first-order transitions 
occur independent of whether the shape of surface is me- 
chanically maintained by the skeleton (= the domain bound- 
ary) or by the surface itself. 

Vertices can freely diffuse inside each compartment on 
the fluid surfaces supported by the skeletons, ft is inter- 
esting to see how fluidity influences the transition of the 
skeleton-supported model. Many interesting problems re- 
main to be studied on the surface model with skeletons. 

This work was supported in part by Grant-in- Aid for 
Scientific Research, No. f5560160 and No. 18560185. 
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